%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Archivo: modelo_llegadas.m
% Genera un histograma, plot q-q y aplica el test chi cuadrado a los
% datos correspondientes al "intervalo entre llegadas al sistema".
clear all;

% Parametros
clases = 7;		% Cant de clases para el histograma
ancho = 0.224/clases;		% Ancho de cada clase
significacion = 0.05;	% Nivel de significacion para el test chi cuadrado

% Cargamos las horas de llegadas.
load llegadasregistro;
n = length(llegadasregistro) - 1;

% Calculamos los intervalos ENTRE llegadas.
for i=1:n
    intervalos(i) = llegadasregistro(i+1) - llegadasregistro(i);
end

% Armamos un histograma
marcas = (0:clases) * ancho + ancho/2;
[f marcas] = hist(intervalos, marcas);

% Como asumimos que los datos provienen de una distribucion 
% exponencial, estimamos lambda como la inversa de la media.
lambda_est = 1 / ((f*marcas') / n)

% Calculamos los valores teoricos de la distribucion, para cada
% una de las marcas de clase.
sup = 1 - exp( -lambda_est * (marcas+ancho/2));
inf = 1 - exp( -lambda_est * (marcas-ancho/2));
teorico = n * (sup - inf);

% Generamos un set de datos con mas valores para poder graficar una curva
% mas suave en el histograma. 
dt = ancho/5;
t = 0:dt:max(marcas);
sup = 1 - exp( -lambda_est * (t+dt/2));
inf = 1 - exp( -lambda_est * (t-dt/2));
z = n * 5 * (sup - inf);

%% histogram for the exponential distribution:
% Mostramos el grafico y lo guardamos en formato eps
%__gnuplot_set__ terminal unknown;

iStdSize = 12;

figure(4);
clf;
hold on;

hist(intervalos, marcas);
plot(t, z, 'k--');

h = findobj(gca,'Type','patch');
set(h,'FaceColor', 'w','EdgeColor','k')
l = legend( '\it Histograma de tiempos entre arribos al sistema', '\it PDF exponencial');  
set( l, 'Interpreter', 'tex', 'Location', 'NorthWest', ...
        'FontName', 'Times', 'FontSize', 9 );
xlbl = xlabel('\it Tiempos entre arribos al sistema [h]');
ylbl = ylabel('\it Frecuancia absoluta');
set(xlbl, 'Interpreter', 'tex', 'FontName', 'Times', 'FontSize', iStdSize);
set(ylbl, 'Interpreter', 'tex', 'FontName', 'Times', 'FontSize', iStdSize);
axis normal;
grid on;

hold off;
%__gnuplot_set__ encoding iso_8859_1;
%__gnuplot_set__ terminal postscript eps;
%__gnuplot_set__ output "../graficos/histograma_llegadas.eps";
%replot;
%closeplot;

%% statistical hypotesis test:
% Calculamos el estadistico chi cuadrado
chi_est = 0;
for i=1:clases
    chi_est = chi_est + ((teorico(i) - f(i))^2 / teorico(i));
end
chi_est

% Calculamos el valor critico 
chi_tabla = chi2inv(1 - significacion, clases - 2)

% Imprimimos resultado del test
if chi_est >= chi_tabla
	fprintf('Los datos no provienen de una distribucion exponencial.\n');
else
	fprintf('No puedo decir que no sean de una distribucion exponencial.\n');
end

% Calculamos quantiles de la muestra y teoricos
quantil_muestra = sort(intervalos);
quantil_teorico = (-1 / lambda_est) * log( 1 - (((1:n) - 0.5) / n));

%% Q-Q Plot of the exponential distribution:
% Los graficamos
%__gnuplot_set__ terminal unknown;

figure(5);
clf;
hold on;

plot(quantil_muestra, quantil_teorico, 'ko', ...
        'LineWidth', 1, 'MarkerEdgeColor', 'k', ...
        'MarkerFaceColor', 'w', 'MarkerSize', 5);
n = 0:dt:0.25;
plot(n, n, '-k', 'LineWidth', 1); 
xlbl = xlabel('\it Cuantiles de la muestra');
ylbl = ylabel('\it Cuantiles del modelo teorico');

set(xlbl, 'Interpreter', 'tex', 'FontName', 'Times', 'FontSize', iStdSize);
set(ylbl, 'Interpreter', 'tex', 'FontName', 'Times', 'FontSize', iStdSize);
axis normal;
grid on;

hold off;
%__gnuplot_set__ encoding iso_8859_1;
%__gnuplot_set__ terminal postscript eps;
%__gnuplot_set__ output "../graficos/plot_qq_llegadas.eps";
%replot;
%closeplot;


%% Aplicamos el test KS:

% Normalizo los datos en el intervalo (0, 1)
llegadasregistro = (llegadasregistro - llegadasregistro(1)) / (llegadasregistro(length(llegadasregistro)) - llegadasregistro(1));
n = length(llegadasregistro);


% Genero un histograma
%gset terminal x11
figure(6);
clf;

hist(llegadasregistro, 10);
h = findobj(gca,'Type','patch');
set(h,'FaceColor', 'w','EdgeColor','k');

[f marcas] = hist(llegadasregistro, clases);

% Computo los estadisticos
f = f / n;
F(1) = 0;
for k =1:clases
	F(k+1) = sum(f(1:k));
end

Dmas = max(F(2:clases+1) - marcas);
Dmenos = max(marcas - F(1:clases));
D = max(Dmas, Dmenos) % max. distance D of the KS-test ...

alpha = 0.05;
% get the critical value ...
c = getCriticalValue_KS(alpha, clases)
% verify the null hypothesis ...
checkNullHypothesis(D, c);

%% results:
% nr. of classes: 7
% 
% lambda_est =
% 
%    18.8253
% 
% 
% chi_est =
% 
%    10.7092
% 
% 
% chi_tabla =
% 
%    11.0705
% 
% No puedo decir que no sean de una distribucion exponencial.
% 
% D =
% 
%     0.1238
% 
% 
% cIdx =
% 
%      4
% 
% 
% n =
% 
%      7
% 
% 
% c =
% 
%     0.4860
% 
% The value does not exceed the critical value c (0.123762 < 0.486000).
% Null hypothesis accepted!
